
setwd("")
library(reshape2)
library(xtable)
library(stargazer)
library(lmtest)
library(DataCombine)
library(plm)
library(Amelia)

data <- read.csv("Beckman_RC_data.csv", header = T)

#View(data)

### Summary statistics (Table 2)
stargazer(data)


### Impute data


set.seed(1234)
data_imp <- amelia(x = data, m=10, ts= "Year", cs="Country",ords=c("checks","execrlc","yrcurnt","yrsoffc","recession_tot","recovery_tot"),
                   noms = c("eu"), polytime = 3,leads = c("sb1"))



############################################
############ MODELS  ######################
############################################

# Mod 1 
b_out_recess <- NULL # Ests
se_out_recess <- NULL # SEs
r_square_recess <- NULL # R2
adj_r_square_recess <-NULL # Adj R2
serial_cor_reces <- NULL # Serial cor Est
serial_cor_p_reces <- NULL # Serial Cor se
# Mod 2
b_out_recov <- NULL
se_out_recov <- NULL
r_square_recov <- NULL
adj_r_square_recov <- NULL
serial_cor_recov <- NULL
serial_cor_p_recov <- NULL
# Mod 3
b_out_recess_dum <- NULL
se_out_recess_dum <- NULL
r_square_recess_dum <- NULL
adj_r_square_reces_dum <- NULL
serial_cor_reces_dum <- NULL
serial_cor_p_reces_dum <- NULL
# Mod 4
b_out_recov_dum <- NULL
se_out_recov_dum <- NULL
r_square_recov_dum <- NULL
adj_r_square_recov_dum <- NULL
serial_cor_recov_dum <- NULL
serial_cor_p_recov_dum <- NULL

# Mod 5
b_out_recess_red <- NULL
se_out_recess_red <- NULL
r_square_recess_red <- NULL
adj_r_square_recess_red <-NULL
serial_cor_reces_red <- NULL
serial_cor_p_reces_red <- NULL
# Mod 6
b_out_recov_red <- NULL
se_out_recov_red <- NULL
r_square_recov_red <- NULL
adj_r_square_recov_red <- NULL
serial_cor_recov_red <- NULL
serial_cor_p_recov_red <- NULL
# Mod 7
b_out_recess_dum_red <- NULL
se_out_recess_dum_red <- NULL
r_square_recess_dum_red <- NULL
adj_r_square_reces_dum_red <- NULL
serial_cor_reces_dum_red <- NULL
serial_cor_p_reces_dum_red <- NULL
# Mod 8
b_out_recov_dum_red <- NULL
se_out_recov_dum_red <- NULL
r_square_recov_dum_red <- NULL
adj_r_square_recov_dum_red <- NULL
serial_cor_recov_dum_red <- NULL
serial_cor_p_recov_dum_red<- NULL



for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  dat$t <- dat$Year - min(as.numeric(paste(dat$Year)))
  dat$t2 <- dat$t^2
  dat$t3 <- dat$t^3
  dat <- pdata.frame(dat, index=c("Country", "Year"))
  dat[,c(4:9,16,19,21:34)] <- scale(dat[,c(4:9,16,19,21:34)]) # Scale vars
  
  #recession_tot model (mod 1)
  femodel_out_reces <- plm(sb1 ~ sbal+sbal_abroad+recession_tot_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                             sbal+sbal_abroad+recession_tot_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  # Standard Errors
  se_recess <-coeftest(femodel_out_reces, vcov=vcovHC(femodel_out_reces,cluster="group"))
  # Get estimates
  b_out_recess <- rbind(b_out_recess, femodel_out_reces$coef)
  se_out_recess <- rbind(se_out_recess, se_recess[,2])
  # Get r2/adj r2
  r_square_recess <- rbind(r_square_recess, summary(femodel_out_reces)$r.square[1])
  adj_r_square_recess <- rbind(adj_r_square_recess, summary(femodel_out_reces)$r.square[2])
  # Get est SE for serial cor
  serial_cor_reces <- rbind(serial_cor_reces, pbgtest(femodel_out_reces,type="Chisq",order=1)$coefficients[22]) # Last est (same with vcov)
  serial_cor_p_reces <- rbind(serial_cor_p_reces, sqrt(pbgtest(femodel_out_reces,type="Chisq",order=1)$vcov[22,22]))
  
  
  # Mod 2
  
  femodel_out_recov <- plm(sb1 ~ sbal+sbal_abroad+recovery_tot_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                             sbal+sbal_abroad+recovery_tot_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recov <-coeftest(femodel_out_recov, vcov=vcovHC(femodel_out_recov,cluster="group"))
  b_out_recov <- rbind(b_out_recov, femodel_out_recov$coef)
  se_out_recov <- rbind(se_out_recov, se_recov[,2])
  r_square_recov <- rbind(r_square_recov, summary(femodel_out_recov)$r.square[1])
  adj_r_square_recov <- rbind(adj_r_square_recov, summary(femodel_out_recov)$r.square[2])
  serial_cor_recov <- rbind(serial_cor_recov, pbgtest(femodel_out_recov,type="Chisq",order=1)$coefficients[22])
  serial_cor_p_recov <- rbind(serial_cor_p_recov, sqrt(pbgtest(femodel_out_recov,type="Chisq",order=1)$vcov[22,22]))
  
  # Mod 3
  femodel_out_reces_dum <- plm(sb1 ~ sbal+sbal_abroad+recession_dum_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                                 sbal+sbal_abroad+recession_dum_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recess_dum <-coeftest(femodel_out_reces_dum, vcov=vcovHC(femodel_out_reces_dum,cluster="group"))
  b_out_recess_dum <- rbind(b_out_recess_dum, femodel_out_reces_dum$coef)
  se_out_recess_dum <- rbind(se_out_recess_dum, se_recess_dum[,2])
  r_square_recess_dum <- rbind(r_square_recess_dum, summary(femodel_out_reces_dum)$r.square[1])
  adj_r_square_reces_dum <- rbind(adj_r_square_reces_dum, summary(femodel_out_reces_dum)$r.square[2])
  serial_cor_reces_dum <- rbind(serial_cor_reces_dum, pbgtest(femodel_out_reces_dum,type="Chisq",order=1)$coefficients[22])
  serial_cor_p_reces_dum <- rbind(serial_cor_p_reces_dum, sqrt(pbgtest(femodel_out_reces_dum,type="Chisq",order=1)$vcov[22,22]))
  
  # Mod 4
  femodel_out_recov_dum <- plm(sb1 ~ sbal+sbal_abroad+recovery_dum_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                                 sbal+sbal_abroad+recovery_dum_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recov_dum <-coeftest(femodel_out_recov_dum, vcov=vcovHC(femodel_out_recov_dum,cluster="group"))
  b_out_recov_dum <- rbind(b_out_recov_dum, femodel_out_recov_dum$coef)
  se_out_recov_dum <- rbind(se_out_recov_dum, se_recov_dum[,2])
  r_square_recov_dum <- rbind(r_square_recov_dum, summary(femodel_out_recov_dum)$r.square[1])
  adj_r_square_recov_dum <- rbind(adj_r_square_recov_dum, summary(femodel_out_recov_dum)$r.square[2])
  serial_cor_recov_dum <- rbind(serial_cor_recov_dum, pbgtest(femodel_out_recov_dum,type="Chisq",order=1)$coefficients[22])
  serial_cor_p_recov_dum <- rbind(serial_cor_p_recov_dum, sqrt(pbgtest(femodel_out_recov_dum,type="Chisq",order=1)$vcov[22,22]))
  
  
  #### REDUCED FORM MODELS
  
  # Mod 5
  femodel_out_reces_red <- plm(sb1 ~ sbal+sbal_abroad+recession_tot_abroad+t+t2+t3|
                                 sbal+sbal_abroad+recession_tot_dist+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  se_recess_red <-coeftest(femodel_out_reces_red, vcov=vcovHC(femodel_out_reces_red,cluster="group"))
  b_out_recess_red <- rbind(b_out_recess_red, femodel_out_reces_red$coef)
  se_out_recess_red <- rbind(se_out_recess_red, se_recess_red[,2])
  r_square_recess_red <- rbind(r_square_recess_red, summary(femodel_out_reces_red)$r.square[1])
  adj_r_square_recess_red <- rbind(adj_r_square_recess_red, summary(femodel_out_reces_red)$r.square[2])
  serial_cor_reces_red <- rbind(serial_cor_reces_red, pbgtest(femodel_out_reces_red,type="Chisq",order=1)$coefficients[7])
  serial_cor_p_reces_red <- rbind(serial_cor_p_reces_red, sqrt(pbgtest(femodel_out_reces_red,type="Chisq",order=1)$vcov[7,7]))
  
  
  # Mod 6
  femodel_out_recov_red <- plm(sb1 ~ sbal+sbal_abroad+recovery_tot_abroad+t+t2+t3|
                                 sbal+sbal_abroad+recovery_tot_dist+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recov_red <-coeftest(femodel_out_recov_red, vcov=vcovHC(femodel_out_recov_red,cluster="group"))
  b_out_recov_red <- rbind(b_out_recov_red, femodel_out_recov_red$coef)
  se_out_recov_red <- rbind(se_out_recov_red, se_recov_red[,2])
  r_square_recov_red <- rbind(r_square_recov_red, summary(femodel_out_recov_red)$r.square[1])
  adj_r_square_recov_red <- rbind(adj_r_square_recov_red, summary(femodel_out_recov_red)$r.square[2])
  serial_cor_recov_red <- rbind(serial_cor_recov_red, pbgtest(femodel_out_recov_red,type="Chisq",order=1)$coefficients[7])
  serial_cor_p_recov_red <- rbind(serial_cor_p_recov_red, sqrt(pbgtest(femodel_out_recov_red,type="Chisq",order=1)$vcov[7,7]))
  
  # Mod 7
  femodel_out_reces_dum_red <- plm(sb1 ~ sbal+sbal_abroad+recession_dum_abroad+t+t2+t3|
                                     sbal+sbal_abroad+recession_dum_dist+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recess_dum_red <-coeftest(femodel_out_reces_dum_red, vcov=vcovHC(femodel_out_reces_dum_red,cluster="group"))
  b_out_recess_dum_red <- rbind(b_out_recess_dum_red, femodel_out_reces_dum_red$coef)
  se_out_recess_dum_red <- rbind(se_out_recess_dum_red, se_recess_dum_red[,2])
  r_square_recess_dum_red <- rbind(r_square_recess_dum_red, summary(femodel_out_reces_dum_red)$r.square[1])
  adj_r_square_reces_dum_red <- rbind(adj_r_square_reces_dum_red, summary(femodel_out_reces_dum_red)$r.square[2])
  serial_cor_reces_dum_red <- rbind(serial_cor_reces_dum_red, pbgtest(femodel_out_reces_dum_red,type="Chisq",order=1)$coefficients[7])
  serial_cor_p_reces_dum_red <- rbind(serial_cor_p_reces_dum_red, sqrt(pbgtest(femodel_out_reces_dum_red,type="Chisq",order=1)$vcov[7,7]))
  
  # Mod 8
  femodel_out_recov_dum_red <- plm(sb1 ~ sbal+sbal_abroad+recovery_dum_abroad+t+t2+t3|
                                     sbal+sbal_abroad+recovery_dum_dist+t+t2+t3, data=dat, model="within", na.rm=TRUE)
  
  se_recov_dum_red <-coeftest(femodel_out_recov_dum_red, vcov=vcovHC(femodel_out_recov_dum_red,cluster="group"))
  b_out_recov_dum_red <- rbind(b_out_recov_dum_red, femodel_out_recov_dum_red$coef)
  se_out_recov_dum_red <- rbind(se_out_recov_dum_red, se_recov_dum_red[,2])
  r_square_recov_dum_red <- rbind(r_square_recov_dum_red, summary(femodel_out_recov_dum_red)$r.square[1])
  adj_r_square_recov_dum_red <- rbind(adj_r_square_recov_dum_red, summary(femodel_out_recov_dum_red)$r.square[2])
  serial_cor_recov_dum_red <- rbind(serial_cor_recov_dum_red, pbgtest(femodel_out_recov_dum_red,type="Chisq",order=1)$coefficients[7])
  serial_cor_p_recov_dum_red <- rbind(serial_cor_p_recov_dum_red, sqrt(pbgtest(femodel_out_recov_dum_red,type="Chisq",order=1)$vcov[7,7]))
  
  
}

# Mod 1
# Combine results from imputations
combined_results1 <- mi.meld(q=b_out_recess, se=se_out_recess)
# Get z scores
mi_z1 <- combined_results1$q.mi/combined_results1$se.mi
# Get p-values
mi_p1 <- 2*(1 - pnorm(abs(mi_z1)))
# Combine the above (minus z-score)
combined_final1 <- cbind(t(combined_results1$q.mi), t(combined_results1$se.mi), t(mi_p1))
#combined_final1 <- cbind(t(round(combined_results1$q.mi,3)), t(round(combined_results1$se.mi,3)), round(t(mi_p1),4))
colnames(combined_final1) <- c("EST", "SE", "P")



# Serial cor estimates for mod 1; do the same as above 
sc_reces <- mi.meld(q=serial_cor_reces, se=serial_cor_p_reces)
mi_zsc_reces <- sc_reces$q.mi/sc_reces$se.mi
mi_psc_reces <- 2*(1 - pnorm(abs(mi_zsc_reces)))
#sc_reces_fin <- t(sc_reces$q.mi)
sc_reces_fin <- cbind(t(sc_reces$q.mi), t(sc_reces$se.mi), t(mi_psc_reces))
#sc_reces_fin <- cbind(t(round(sc_reces$q.mi,3)), t(round(sc_reces$se.mi,3)), round(t(mi_psc_reces),4))
colnames(sc_reces_fin) <- c("EST", "SE", "P")
#sc_reces_fin


# Mod 2
combined_results2 <- mi.meld(q=b_out_recov, se=se_out_recov)
mi_z2 <- combined_results2$q.mi/combined_results2$se.mi
mi_p2 <- 2*(1 - pnorm(abs(mi_z2)))
combined_final2 <- cbind(t(combined_results2$q.mi), t(combined_results2$se.mi), t(mi_p2))
#combined_final2 <- cbind(t(round(combined_results2$q.mi,3)), t(round(combined_results2$se.mi,3)), round(t(mi_p2),4))
colnames(combined_final2) <- c("EST", "SE", "P")

sc_recov <- mi.meld(q=serial_cor_recov, se=serial_cor_p_recov)
mi_zsc_recov <- sc_recov$q.mi/sc_recov$se.mi
mi_psc_recov <- 2*(1 - pnorm(abs(mi_zsc_recov)))
sc_recov_fin <- cbind(t(sc_recov$q.mi), t(sc_recov$se.mi), t(mi_psc_recov))
#sc_recov_fin <- cbind(t(round(sc_recov$q.mi,3)), t(round(sc_recov$se.mi,3)), round(t(mi_psc_recov),4))
colnames(sc_recov_fin) <- c("EST", "SE", "P")
#combined_final2
#combined_final2 <- rbind(combined_final2, sc_recov_fin)


# Mod 3

combined_results3 <- mi.meld(q=b_out_recess_dum, se=se_out_recess_dum)
mi_z3 <- combined_results3$q.mi/combined_results3$se.mi
mi_p3 <- 2*(1 - pnorm(abs(mi_z3)))
combined_final3 <- cbind(t(combined_results3$q.mi), t(combined_results3$se.mi), t(mi_p3))
#combined_final3 <- cbind(t(round(combined_results3$q.mi,3)), t(round(combined_results3$se.mi,3)), round(t(mi_p3),4))
colnames(combined_final3) <- c("EST", "SE", "P")

sc_reces_dum <- mi.meld(q=serial_cor_reces_dum, se=serial_cor_p_reces_dum)
mi_zsc_reces_dum <- sc_reces_dum$q.mi/sc_reces$se.mi
mi_psc_reces_dum <- 2*(1 - pnorm(abs(mi_zsc_reces_dum)))
sc_reces_dum_fin <- cbind(t(sc_reces_dum$q.mi), t(sc_reces_dum$se.mi), t(mi_psc_reces_dum))
#sc_reces_dum_fin <- cbind(t(round(sc_reces_dum$q.mi,3)), t(round(sc_reces_dum$se.mi,3)), round(t(mi_psc_reces_dum),4))
colnames(sc_reces_dum_fin) <- c("EST", "SE", "P")

#combined_final3 <- rbind(combined_final3, sc_reces_dum_fin)

# Mod 4

combined_results4 <- mi.meld(q=b_out_recov_dum, se=se_out_recov_dum)
mi_z4 <- combined_results4$q.mi/combined_results4$se.mi
mi_p4 <- 2*(1 - pnorm(abs(mi_z4)))
combined_final4 <- cbind(t(combined_results4$q.mi), t(combined_results4$se.mi), t(mi_p4))
#combined_final4 <- cbind(t(round(combined_results4$q.mi,3)), t(round(combined_results4$se.mi,3)), round(t(mi_p4),4))
colnames(combined_final4) <- c("EST", "SE", "P")

sc_recov_dum <- mi.meld(q=serial_cor_recov_dum, se=serial_cor_p_recov_dum)
mi_zsc_recov_dum <- sc_recov_dum$q.mi/sc_recov_dum$se.mi
mi_psc_recov_dum <- 2*(1 - pnorm(abs(mi_zsc_recov_dum)))
sc_recov_dum_fin <- cbind(t(sc_recov_dum$q.mi), t(sc_recov_dum$se.mi), t(mi_psc_recov_dum))
#sc_recov_dum_fin <- cbind(t(round(sc_recov_dum$q.mi,3)), t(round(sc_recov_dum$se.mi,3)), round(t(mi_psc_recov_dum),4))
colnames(sc_recov_dum_fin) <- c("EST", "SE", "P")

#combined_final4 <- rbind(combined_final4, sc_recov_dum_fin)

# Mod 5
combined_results5 <- mi.meld(q=b_out_recess_red, se=se_out_recess_red)
mi_z5 <- combined_results5$q.mi/combined_results5$se.mi
mi_p5 <- 2*(1 - pnorm(abs(mi_z5)))
combined_final5 <- cbind(t(combined_results5$q.mi), t(combined_results5$se.mi), t(mi_p5))
#combined_final5 <- cbind(t(round(combined_results5$q.mi,3)), t(round(combined_results5$se.mi,3)), round(t(mi_p5),4))
colnames(combined_final5) <- c("EST", "SE", "P")


sc_reces_red <- mi.meld(q=serial_cor_reces_red, se=serial_cor_p_reces_red)
mi_zsc_reces_red <- sc_reces_red$q.mi/sc_reces_red$se.mi
mi_psc_reces_red <- 2*(1 - pnorm(abs(mi_zsc_reces_red)))
sc_reces_fin_red <- cbind(t(sc_reces_red$q.mi), t(sc_reces_red$se.mi), t(mi_psc_reces_red))
#sc_reces_fin_red <- cbind(t(round(sc_reces_red$q.mi,3)), t(round(sc_reces_red$se.mi,3)), round(t(mi_psc_reces_red),4))
colnames(sc_reces_fin_red) <- c("EST", "SE", "P")
#sc_reces_fin_red


#combined_final5 <- rbind(combined_final5, sc_reces_fin_red)

# Mod 6
combined_results6 <- mi.meld(q=b_out_recov_red, se=se_out_recov_red)
mi_z6 <- combined_results6$q.mi/combined_results6$se.mi
mi_p6 <- 2*(1 - pnorm(abs(mi_z6)))
combined_final6 <- cbind(t(combined_results6$q.mi), t(combined_results6$se.mi), t(mi_p6))
#combined_final6 <- cbind(t(round(combined_results6$q.mi,3)), t(round(combined_results6$se.mi,3)), round(t(mi_p6),4))
colnames(combined_final6) <- c("EST", "SE", "P")
#combined_final6


sc_recov_red <- mi.meld(q=serial_cor_recov_red, se=serial_cor_p_recov_red)
mi_zsc_recov_red <- sc_recov_red$q.mi/sc_recov_red$se.mi
mi_psc_recov_red <- 2*(1 - pnorm(abs(mi_zsc_recov_red)))
sc_recov_fin_red <- cbind(t(sc_recov_red$q.mi), t(sc_recov_red$se.mi), t(mi_psc_recov_red))
#sc_recov_fin_red <- cbind(t(round(sc_recov_red$q.mi,3)), t(round(sc_recov_red$se.mi,3)), round(t(mi_psc_recov_red),4))
colnames(sc_recov_fin_red) <- c("EST", "SE", "P")
#sc_recov_fin_red

#combined_final6 <- rbind(combined_final6, sc_recov_fin_red)


# Mod 7

combined_results7 <- mi.meld(q=b_out_recess_dum_red, se=se_out_recess_dum_red)
mi_z7 <- combined_results7$q.mi/combined_results7$se.mi
mi_p7 <- 2*(1 - pnorm(abs(mi_z7)))
combined_final7 <- cbind(t(combined_results7$q.mi), t(combined_results7$se.mi), t(mi_p7))
#combined_final7 <- cbind(t(round(combined_results7$q.mi,3)), t(round(combined_results7$se.mi,3)), round(t(mi_p7),4))
colnames(combined_final7) <- c("EST", "SE", "P")



sc_reces_dum_red <- mi.meld(q=serial_cor_reces_dum_red, se=serial_cor_p_reces_dum_red)
mi_zsc_reces_dum_red <- sc_reces_dum_red$q.mi/sc_reces_dum_red$se.mi
mi_psc_reces_dum_red <- 2*(1 - pnorm(abs(mi_zsc_reces_dum_red)))
sc_reces_dum_fin_red <- cbind(t(sc_reces_dum_red$q.mi), t(sc_reces_dum_red$se.mi), t(mi_psc_reces_dum))
#sc_reces_dum_fin_red <- cbind(t(round(sc_reces_dum_red$q.mi,3)), t(round(sc_reces_dum_red$se.mi,3)), round(t(mi_psc_reces_dum),4))
colnames(sc_reces_dum_fin_red) <- c("EST", "SE", "P")

#combined_final7 <- rbind(combined_final7, sc_reces_dum_fin_red)


# Mod 8

combined_results8 <- mi.meld(q=b_out_recov_dum_red, se=se_out_recov_dum_red)
mi_z8 <- combined_results8$q.mi/combined_results8$se.mi
mi_p8 <- 2*(1 - pnorm(abs(mi_z8)))
combined_final8 <- cbind(t(combined_results8$q.mi), t(combined_results8$se.mi), t(mi_p8))
#combined_final8 <- cbind(t(round(combined_results8$q.mi,3)), t(round(combined_results8$se.mi,3)), round(t(mi_p8),4))
colnames(combined_final8) <- c("EST", "SE", "P")

sc_recov_dum_red <- mi.meld(q=serial_cor_recov_dum_red, se=serial_cor_p_recov_dum_red)
mi_zsc_recov_dum_red <- sc_recov_dum_red$q.mi/sc_recov_dum_red$se.mi
mi_psc_recov_dum_red <- 2*(1 - pnorm(abs(mi_zsc_recov_dum_red)))
sc_recov_dum_fin_red <- cbind(t(sc_recov_dum_red$q.mi), t(sc_recov_dum_red$se.mi), t(mi_psc_recov_dum_red))
#sc_recov_dum_fin_red <- cbind(t(round(sc_recov_dum_red$q.mi,3)), t(round(sc_recov_dum_red$se.mi,3)), round(t(mi_psc_recov_dum_red),4))
colnames(sc_recov_dum_fin_red) <- c("EST", "SE", "P")

#combined_final8 <- rbind(combined_final8, sc_recov_dum_fin_red)


## Get r-square, adjusted r-square

r2_1 <- sum(r_square_recess[,1])/data_imp$m
r2_adj_1 <- sum(adj_r_square_recess[,1])/data_imp$m
r2_2 <- sum(r_square_recov[,1])/data_imp$m
r2_adj_2 <- sum(adj_r_square_recov[,1])/data_imp$m
r2_3 <- sum(r_square_recess_dum[,1])/data_imp$m
r2_adj_3 <- sum(adj_r_square_reces_dum[,1])/data_imp$m
r2_4 <- sum(r_square_recov_dum[,1])/data_imp$m
r2_adj_4 <- sum(adj_r_square_recov_dum[,1])/data_imp$m
r2_5 <- sum(r_square_recess_red[,1])/data_imp$m
r2_adj_5 <- sum(adj_r_square_recess_red[,1])/data_imp$m
r2_6 <- sum(r_square_recov_red[,1])/data_imp$m
r2_adj_6 <- sum(adj_r_square_recov_red[,1])/data_imp$m
r2_7 <- sum(r_square_recess_dum_red[,1])/data_imp$m
r2_adj_7 <- sum(adj_r_square_reces_dum_red[,1])/data_imp$m
r2_8 <- sum(r_square_recov_dum_red[,1])/data_imp$m
r2_adj_8 <- sum(adj_r_square_recov_dum_red[,1])/data_imp$m

#########################################################
#################### Now Print Models ######
#########################################################

coefs <- list(combined_final1[,1],combined_final2[,1],combined_final3[,1],combined_final4[,1],
              combined_final5[,1],combined_final6[,1],combined_final7[,1],combined_final8[,1])

ses <- list(combined_final1[,2],combined_final2[,2],combined_final3[,2],combined_final4[,2],
            combined_final5[,2],combined_final6[,2],combined_final7[,2],combined_final8[,2])

pvals <- list(combined_final1[,3],combined_final2[,3],combined_final3[,3],combined_final4[,3],
              combined_final5[,3],combined_final6[,3],combined_final7[,3],combined_final8[,3])


stargazer(femodel_out_reces,femodel_out_recov,femodel_out_reces_dum,femodel_out_recov_dum,
          femodel_out_reces_red,femodel_out_recov_red,femodel_out_reces_dum_red,
          femodel_out_recov_dum_red, 
          coef = coefs, se= ses, p= pvals,
          omit.stat=c("f"), 
          dep.var.labels = c("Adjusted Budget Balance$_{t+1}$"),
          covariate.labels=c("Primary Balance","Spatial Primary Balance","Recession Abroad (Tot. Quarters)","Recovery Abroad (Tot. Quarters)",
                             "Recession Abroad (Prop.)","Recovery Abroad (Prop.)","Change GDP","Domestic Recession", "Domestic Recovery",
                             "Trade Openness","Inflation","GDP Per Capita (Log)","China Growth","Interest Rates","Debt to GDP","EU","EU times Debt",
                             "Vetoes","Years Exec. in Office","Years to Next Exec. Elec.","Legislative Frac."),
          notes = "Cluster robust standard errors in parentheses (clustered by country)"
)

# Get AR(1) terms

scs <- list(sc_reces_fin[,1], sc_recov_fin[,1],sc_reces_dum_fin[,1],sc_recov_dum_fin[,1],
             sc_reces_fin_red[,1],sc_recov_fin_red[,1], sc_reces_dum_fin_red[,1],sc_recov_dum_fin_red[,1])
scs_se <- list(sc_reces_fin[,2], sc_recov_fin[,2],sc_reces_dum_fin[,2],sc_recov_dum_fin[,2],
                sc_reces_fin_red[,2],sc_recov_fin_red[,2], sc_reces_dum_fin_red[,2],sc_recov_dum_fin_red[,2])

scs_p <- list(sc_reces_fin[,3], sc_recov_fin[,3],sc_reces_dum_fin[,3],sc_recov_dum_fin[,3],
               sc_reces_fin_red[,3],sc_recov_fin_red[,3], sc_reces_dum_fin_red[,3],sc_recov_dum_fin_red[,2])

# need placeholder model 
mod_1 <- lm(sbal~-1+ recession_tot_abroad, data=data)

stargazer(mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,
          coef=scs, se=scs_se, p = scs_p, omit.stat = 'F',
          covariate.labels = 'BG AR(1) $Phi$')

# Just take estimate and ses from above!




r2 <- cbind(r2_1,r2_2,r2_3,r2_4,r2_5,r2_6,r2_7,r2_8)
adj_r2 <- cbind(r2_adj_1,r2_adj_2,r2_adj_3,r2_adj_4,r2_adj_5,r2_adj_6,r2_adj_7,r2_adj_8)

xtable(r2)
xtable(adj_r2)











##### LISTWISE DELETION

data <- read.csv("Beckman_RC_data.csv", header = T)

data$t <- data$Year - min(as.numeric(paste(data$Year)))
data$t2 <- data$t^2
data$t3 <- data$t^3

data <- pdata.frame(data, index=c("Country", "Year"))

femodel_out_reces <- plm(sb1 ~ sbal+sbal_abroad+recession_tot_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                           sbal+sbal_abroad+recession_tot_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=data, model="within", na.rm=TRUE)

# Standard Errors
se1 <- sqrt(diag(vcovHC(femodel_out_reces,cluster="group")))
# Serial cor estimate
sc_est1 <- pbgtest(femodel_out_reces,type="Chisq",order=1)$coefficients[22][[1]] # Last est (same with vcov)
# Serial cor standard error
sc_se1 <- sqrt(pbgtest(femodel_out_reces,type="Chisq",order=1)$vcov[22,22])
# Serial cor p
sc_p1 <- pbgtest(femodel_out_recov,type="Chisq",order=1)[4][[1]]


# Mod 2

femodel_out_recov <- plm(sb1 ~ sbal+sbal_abroad+recovery_tot_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                           sbal+sbal_abroad+recovery_tot_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=data, model="within", na.rm=TRUE)

se2 <- sqrt(diag(vcovHC(femodel_out_recov,cluster="group")))

sc_est2 <- pbgtest(femodel_out_recov,type="Chisq",order=1)$coefficients[22][[1]]
sc_se2 <- sqrt(pbgtest(femodel_out_recov,type="Chisq",order=1)$vcov[22,22])
sc_p2 <- pbgtest(femodel_out_recov,type="Chisq",order=1)[4][[1]]


# Mod 3
femodel_out_reces_dum <- plm(sb1 ~ sbal+sbal_abroad+recession_dum_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                               sbal+sbal_abroad+recession_dum_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=data, model="within", na.rm=TRUE)

se3 <- sqrt(diag(vcovHC(femodel_out_reces_dum,cluster="group")))

sc_est3 <- pbgtest(femodel_out_reces_dum,type="Chisq",order=1)$coefficients[22][[1]]
sc_se3 <- sqrt(pbgtest(femodel_out_reces_dum,type="Chisq",order=1)$vcov[22,22])
sc_p3 <- pbgtest(femodel_out_reces_dum,type="Chisq",order=1)[4][[1]]


# Mod 4
femodel_out_recov_dum <- plm(sb1 ~ sbal+sbal_abroad+recovery_dum_abroad+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3|
                               sbal+sbal_abroad+recovery_dum_dist+change_GDP+recession_tot+recovery_tot+openness+inflation+log_gdp+China_pct+int_spread+gross_debt+eu+euro_debt+checks+yrcurnt+yrsoffc+maj+t+t2+t3, data=data, model="within", na.rm=TRUE)

se4 <- sqrt(diag(vcovHC(femodel_out_recov_dum,cluster="group")))

sc_est4 <- pbgtest(femodel_out_recov_dum,type="Chisq",order=1)$coefficients[22][[1]]
sc_se4 <- sqrt(pbgtest(femodel_out_recov_dum,type="Chisq",order=1)$vcov[22,22])
sc_p4 <- pbgtest(femodel_out_recov_dum,type="Chisq",order=1)[4][[1]]


#### REDUCED FORM MODELS

# Mod 5
femodel_out_reces_red <- plm(sb1 ~ sbal+sbal_abroad+recession_tot_abroad+t+t2+t3|
                               sbal+sbal_abroad+recession_tot_dist+t+t2+t3, data=data, model="within", na.rm=TRUE)
se5 <- sqrt(diag(vcovHC(femodel_out_reces_red,cluster="group")))

sc_est5 <- pbgtest(femodel_out_reces_red,type="Chisq",order=1)$coefficients[7][[1]]
sc_se5 <- sqrt(pbgtest(femodel_out_reces_red,type="Chisq",order=1)$vcov[7,7])
sc_p5 <- pbgtest(femodel_out_reces_red,type="Chisq",order=1)[4][[1]]


# Mod 6
femodel_out_recov_red <- plm(sb1 ~ sbal+sbal_abroad+recovery_tot_abroad+t+t2+t3|
                               sbal+sbal_abroad+recovery_tot_dist+t+t2+t3, data=data, model="within", na.rm=TRUE)

se6 <- sqrt(diag(vcovHC(femodel_out_recov_red,cluster="group")))

sc_est6 <- pbgtest(femodel_out_recov_red,type="Chisq",order=1)$coefficients[7][[1]]
sc_se6 <- sqrt(pbgtest(femodel_out_recov_red,type="Chisq",order=1)$vcov[7,7])
sc_p6 <- pbgtest(femodel_out_recov_red,type="Chisq",order=1)[4][[1]]


# Mod 7
femodel_out_reces_dum_red <- plm(sb1 ~ sbal+sbal_abroad+recession_dum_abroad+t+t2+t3|
                                   sbal+sbal_abroad+recession_dum_dist+t+t2+t3, data=data, model="within", na.rm=TRUE)

se7 <- sqrt(diag(vcovHC(femodel_out_reces_dum_red,cluster="group")))

sc_est7 <- pbgtest(femodel_out_reces_dum_red,type="Chisq",order=1)$coefficients[7][[1]]
sc_se7 <- sqrt(pbgtest(femodel_out_reces_dum_red,type="Chisq",order=1)$vcov[7,7])
sc_p7 <- pbgtest(femodel_out_reces_dum_red,type="Chisq",order=1)[4][[1]]


# Mod 8
femodel_out_recov_dum_red <- plm(sb1 ~ sbal+sbal_abroad+recovery_dum_abroad+t+t2+t3|
                                   sbal+sbal_abroad+recovery_dum_dist+t+t2+t3, data=data, model="within", na.rm=TRUE)

se8 <- sqrt(diag(vcovHC(femodel_out_recov_dum_red,cluster="group")))

sc_est8 <- pbgtest(femodel_out_recov_dum_red,type="Chisq",order=1)$coefficients[7][[1]]
sc_se8 <- sqrt(pbgtest(femodel_out_recov_dum_red,type="Chisq",order=1)$vcov[7,7])
sc_p8 <- pbgtest(femodel_out_recov_dum_red,type="Chisq",order=1)[4][[1]]


ses2 <- list(se1,se2,se3,se4,se5,se6,se7,se8)



# Print Main model
stargazer(femodel_out_reces,femodel_out_recov,femodel_out_reces_dum,femodel_out_recov_dum,
          femodel_out_reces_red,femodel_out_recov_red,femodel_out_reces_dum_red,
          femodel_out_recov_dum_red, 
           se= ses2,
          omit.stat=c("f"), 
          dep.var.labels = c("Adjusted Budget Balance$_{t+1}$"),
          covariate.labels=c("Primary Balance","Spatial Primary Balance","Recession Abroad (Tot. Quarters)","Recovery Abroad (Tot. Quarters)",
                             "Recession Abroad (Prop.)","Recovery Abroad (Prop.)","Change GDP","Domestic Recession", "Domestic Recovery",
                             "Trade Openness","Inflation","GDP Per Capita (Log)","China Growth","Interest Rates","Debt to GDP","EU","EU times Debt",
                             "Vetoes","Years Exec. in Office","Years to Next Exec. Elec.","Legislative Frac."),
          notes = "Cluster robust standard errors in parentheses (clustered by country)"
)



# AR(1) ests/ses

sc_ests <- list(sc_est1,sc_est2,sc_est3,sc_est4,sc_est5,sc_est6,sc_est7,sc_est8)
sc_ses <- list(sc_se1,sc_se2,sc_se3,sc_se4,sc_se5,sc_se6,sc_se7,sc_se8)
sc_p <- list(sc_p1,sc_p2,sc_p3,sc_p4,sc_p5,sc_p6,sc_p7,sc_p8)

# Placeholder for table below
mod_1 <- lm(sbal~-1+ recession_tot_abroad, data=data)

# Only take the estimates and ses from table below
stargazer(mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,mod_1,
          coef=sc_ests, se=sc_ses,p=sc_p, omit.stat = 'F',
          covariate.labels = 'BG AR(1) $Phi$')


